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Abstract. ML(n)BiCGStab is a Krylov subspace method for the solution of large, sparse and 
non-symmetric linear systems. In theory, it is a method that lies between the well-known 
BiCGStab and GMRES/FOM. In fact, when n = 1, ML(l)BiCGStab is BiCGStab and when 
n = N, ML(iV)BiCGStab is GMRES/FOM where N is the size of the linear system. There- 
fore, ML(n)BiCGStab is a bridge that connects the Lanczos-based BiCGStab and the Arnoldi- 
based GMRES/FOM. In computation, ML(n)BiCGStab can be much more stable and con- 
verge much faster than BiCGStab when a problem with ill-condition is solved. We have tested 
ML(n)BiCGStab on the standard oil reservoir simulation test data called SPE9 and found that 
ML(n)BiCGStab reduced the total computational time by more than 60% when compared to 
BiCGStab. Tests made on the data from Matrix Market also support the superiority of ML(n)Bi- 
CGStab over BiCGStab. Because of the 0(N 2 ) storage requirement in the full GMRES, one has 
pH _ to adopt a restart strategy to get the storage under control when GMRES is implemented. In 

comparison, ML(n)BiCGStab is a method with only 0(nN) storage requirement and therefore it 
does not need a restart strategy. In this paper, we introduce ML(n)BiCGStab (in particular, a new 
algorithm involving A-transpose), its relations to some existing methods and its implementations. 
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1. Introduction 

ML(n)BiCGStab is a Krylov subspace method for the solution of the linear system 

VO ■ (1.1) Ax = b 

where A € C NxN and b G C . It was introduced by Yeung and Chan [12j in 1999 and its 
algorithms were recently reformulated by Yeung [11]. ML(n)BiCGStab is a natural generalization 
of BiCGStab by van der Vorst [8J, built on the multiple starting Lanczos process rather than 
on the single starting Lanczos process. Its derivation relies on the techniques introduced by 
Sonneveld [6] and van der Vorst [8] in the construction of CGS and BiCGStab. There have 
been three algorithms associated with the ML(n)BiCGStab method so far, depending on how the 
residual vector is defined and whether or not the Hermitian transpose A^ is used. In this 
paper, we shall simply introduce the algorithms and address some implementation issues. For 
more detailed, one is referred to |llj . 

Other extensions of BiCGStab exist. Among them are BiCGStab2 by Gutknecht [9J , BiCGStab(Z) 
by Sleijpen and Fokkema [4] and CPBi-CG by Zhang |13j . 

The outline of the paper is as follows. In $2j we introduce index functions which are help- 
ful in presenting the ML(n)BiCGStab algorithms. In [j3j we present the ML(n)BiCG algo- 
rithm from |12j . from which ML(n)BiCGStab algorithms were derived. In $51 we introduce the 
ML(n)BiCGStab algorithms and their relationships with some existing methods. In imple- 
mentation issues are addressed and conclusions are made in 

2. Index Functions 

Let be given a positive integer n. For all integers k, we define 

9n(k) = [{k - l)/nj and r n (k) = k - ng n {k) 

l 
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where \_ ■ J rounds its argument to the nearest integer towards minus infinity. We call g n and r n 
index functions; they are defined on Z, the set of all integers, with ranges Z and {1,2, • • • ,n}, 
respectively. 

Table 2.1. Simple illustration of the index functions for n = 3. 



k 





12 3 


4 5 6 


7 8 9 


10 11 12 ••• 


9n(k) 


-1 





1 1 1 


2 2 2 


3 3 3 


r n (k) 


3 


1 2 3 


1 2 3 


12 3 


12 3 



If we write 

(2.1) k = jn + i 

with 1 < i < n and j € Z, then 

9n{jn + i)= j and r n {jn + i) = i. 
Table [27T1 illustrates the behavior of g n and r n with n = 3. 

3. A ML(n)BiCG Algorithm 

Parallel to the derivation of BiCGStab from BiCG by Fletcher [I], ML(n)BiCGStab was derived 
from a BiCG-like method named ML(n)BiCG, which was constructed based on the multiple 
starting Lanczos process with n left starting vectors and a single right starting vector. 

Let be given n vectors qi, . . . ,q n E C , which we call left starting vectors or shadow vectors. 
Set 

(3-1) Vk = (A H ) 9 - {k \ rn{k) , k = 1,2,3,... . 

The following algorithm for the solution of eqn (jl.ip is from |12| . 

Algorithm 3.1. ML(n)BiCG 

1. Choose an initial guess xo and n vectors qi,q2, ■ ■ ■ ,qn- 

2. Compute r o = b — Axo and set Pi = qi, go = ?o. 

3. For k = 1, 2, 3, • • • , until convergence: 
4- u k = pf r fe _i/p^Ag fc _i; 

5. x fc =x fc _i + Qtfegfc-i; 

6- r k = r fc _i - a fc Ag fc _i; 

7. For s = max(/c — n, 0), • • • , k — 1 

*• ^ = -P?+1 A + n^ax( fe -n,0) ^ft) /pf + lAg s ; 

9. End 

10. g k = r k + Xy s =max(fc-n,0) ^ s Ss/ 

11. Compute Pfc+i according to eqn $3.1\) 

12. End 

Even though the algorithm has not been tested, it is believed to be numerically instable because 
of Line 11 in which the shadow vectors are repeatedly multiplied by A^, a type of operation which 
is highly sensitive to round-off errors. The algorithm has been introduced only for the purpose of 
developing ML(n)BiCGStab algorithms. 

Relations to some other methods: 



AN INTRODUCTION TO ML(n)BICGSTAB 



3 



Table 4.1. Average cost per iteration of the first ML(ra)BiCGStab algorithm and 
its storage. 



Preconditioning M _1 v 




u ± v, av 


5 

max(4 , 0) 


Matvec Av 


n. 


Saxpy u + av 


22 I 

max(2.5ra + 0.5 + -,6) 
n 


dot product u^v 


2 

n+ 1 + - 

n 


Storage 


A + M+ 

(4n + A)N + 0(n) 



(1) Relation with FOM by Saad and Schultz J3J/. Consider the case where n > N. If we choose 
Qfc = ?fc-i in Algorithm 13 . 1 1 (it is possible since r k -i is computed before is used in Line 
11), then Algorithm 13. II is a FOM algorithm. 

(2) Relation with GMRES by Saad and Schultz |5j/. Consider the case where n > N . If we 
choose qfc = Ar^i in Algorithm 13 .1\ then Algorithm 13.11 is a GMRES algorithm. 

(3) Relation with BiCG. When n = 1, Algorithm 13.11 is a BiCG algorithm. 

4. ML(n)BiCGSTAB Algorithms 

There are three algorithms for the ML(n)BiCGStab method. All were derived from Algorithm 
13.11 The first two algorithms do not involve in their implementation and can be found in |11| . 
The third one, however, needs A^ and is new. Therefore, we spend more space here on the the 
third algorithm. 



4.1. First Algorithm. Let f2&(A) be the polynomial of degree k defined by 

1 if k = 

(1 - w fc A)f4_i(A) if > 0. 

If we define the ML(n)BiCGStab residual r^ by 

SV(fc)+i( A )rfc, tf&>l, 
r , if k = 0, 



rfe 



then Algorithm 13.11 will lead to the first ML(n)BiCGStab algorithm (Algorithm 4.1 in |llj). 
Computational and storage cost based on its preconditioned version (Algorithm 9.1 in |llj ) is 
presented in Table I4TT1 

Relations to some other methods: this first algorithm is a BiCGStab algorithm when n = 1. 
4.2. Second Algorithm. If we define the ML(n)BiCGStab residual r^ by 

l4 - ij k ~\r , if fe = 0, 

then Algorithm 13.11 will lead to the second ML(n)BiCGStab algorithm (Algorithm 5.1 in [llj). 
Computational and storage cost based on its preconditioned version (Algorithm 9.2 in |11| ) is 
presented in Table 1131 

Relations to some other methods: 

(1) Relation with FOM. Consider the case where n > N. If we choose q& = r^_i, then this 
algorithm is a FOM algorithm. 

(2) Relation with GMRES. Consider the case where n > N. If we choose q^ = Ar^_i, then 
this algorithm is a GMRES algorithm. 

(3) Relation with BiCGStab. When n = 1, this algorithm is a BiCGStab algorithm. 
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Table 4.2. Average cost per iteration of the second ML(n)BiCGStab algorithm 
and its storage. 



Preconditioning M _1 v 


i + A 


u ± v, av 


1 


Matvec Av 


n 


Saxpy u + av 


2n + 2 + - 

n 


dot product ir^v 


1 2 
n + 1 + - 

n 


Storage 


A + M+ 

(3n + 5)N + 0(n) 



(4) Relation with IDR(s) by Sonneveld and van Gijzen [7|l 10$. This algorithm is a IDR(n) 
algorithm. 

4.3. Third Algorithm. If we define the ML(n)BiCGStab residual r^ by eqn (|4.ip and get A^ 
involved in its implementation, then through the derivation stages #5 - #8 in Algorithm 13. II 
will lead to the following ML(n)BiCGStab algorithm which we name ML(n)BiCGStabt, standing 
for ML(n)BiCGStab with A-transpose. 



Algorithm 4.1. ML(n)BiCGStabt without preconditioning 

1. Choose an initial guess xo and n vectors qi,q2, ■ ■ ■ , qn- 

2. Compute [fi , ■ ■ ■ , f n -i] = A H [q x , • • • , q n _i] . 

3. Compute vq = b — Axo and go = ro, wo = Ago, Co = qf^wo- 
4- For k = 1, 2, • • • , until convergence: 



5. 


a k 






6. 


If 


r n {k) < n 




7. 
8. 
9. 




Xfc = x fc __i + a k gk-i; 
z w = r k ; gfc = 0; 
For s = max(fc — n, 0) 


rfc = 
' ' ' ) 


10. 






/c s ; 


11. 




z w =z w + 'w s ; 




12. 








13. 




End 




14. 




Sk ~ Zw oj gn{k+1) Ek, 




15. 




For s = g n {k)n, ■■■ ,k 


- 1 


16. 








17. 




gfc = gfc + p s g s ; 




18. 




End 




19. 


Else 




20. 




xfc = x fc _i + a fc gfc_i; 




21. 

22. 




Ufc = rjfc_i - a fc Wfc_i; 
^„(fc+i) = (Au k ) H u k / 


||Au 


23. 




Xfc = x fc + u; fln (fc +1 )Ufc; 


Tfc = 


21 




%w = rfc; gfc = 0; 




25. 




For s = g n (k)n,--- , k 


- 1 


26. 






/c s ; 


27. 




z w = z w + (ii ; w s ; 





rfc-i - afcWfc_i 



(k)n - 1 



■k 1 1 1 t 



_w <?„(fc+l) Au fc + u fc/ 

% fi^ = —U}g„(k+l)Ps k 
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Table 4.3. Average cost per iteration of preconditioned ML(n)BiCGStabt and 
its storage. This table does not count the cost in Lines 1-2 of Algorithm 17.11 



Preconditioning M _1 v 


i + i 


u ± v, av 


1 


Matvec Av 


n 


Saxpy u + av 


1.5n + 2.5 + - 

n 


dot product ir^v 


n+ 1 + - 

n 


Storage 


A + M+ 

(4n + 4)AT + 0(n) 



28. g fc = 

29. End 

30. g fc = z w 

31. End 

32. w fc = Ag fc ; 

33. End 



Cfc = q^ (fc+1) Wfc; 



A preconditioned version of Algorithm 14.11 can be obtained by applying it to AM _1 y = b, then 
recovering x through x = M _1 y. The resulting preconditioned algorithm and its Matlab code 
are attached in JO Computational and storage cost is presented in Table 14.31 

Relations to some other methods: Algorithm 14. II is a BiCGStab algorithm when n = 1. 

5. Implementation Issues 

The following test data were downloaded from Matrix Market. More experiments can be found 
in P3K2]. 

(1) utm5940, TOKAMAK Nuclear Physics (Plasmas). utm5940 contains a 5940 x 5940 real 
unsymmetric matrix A with 83, 842 nonzero entries and a real right-hand side b. 

(2) qc2534, H2PLUS Quantum Chemistry, NEP Collection. qc2534 contains a 2534 x 2534 
complex symmetric indefinite matrix with 463, 360 nonzero entries, but does not provide 
the right-hand side b. We set b = Al with 1 = [1, , 1, • • • , 1] T . 

All computing in this section was done in Matlab Version 7.1 on a Windows XP machine with 
a Pentium 4 processor. ILU(0) preconditioners (p. 294, [2]) were used, initial guess was xo = 
and the stopping criterion was 

||rfc||2/!|b|| 2 < HT 7 

where was the computed residual. Shadow vectors Q = [qi,q2,-"" >Qn] were chosen to be 
Q = [ro, randn(N, n — 1)] for utm5940 and Q = [ro, randn(N, n — 1) +sqrt{— 1) *randn(N, n — 1)] 
for qc2534. 

For the convenience of our presentation, let us introduce the following functions: 

(a) T conv (n) is the time that a ML(n)BiCGStab algorithm takes to converge. 

(b) E(n) = ||b — Ax||2/||b||2 is the true relative error of x where x is the computed solution 
output by a ML(n)BiCGStab algorithm when it converges. 

5.1. Stability. The graphs of E(n) are plotted in Figure ISTTl It can be seen that the computed 
rfc by the second algorithm can easily diverges from its exact counterpart b — Axfc. This diver- 
gence becomes significant when n > 4 for utm5940. By contrast, the computed relative errors 
|j r fc||2/||b||2 by the first and the third algorithms well approximate their corresponding true ones. 
Thus, from this point of view, we consider that the first and the third algorithms are numerically 
more stable than the second algorithm. 
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Figure 5.1. Graphs of E(n) against n. First algorithm: x-mark; Second algo- 
rithm: o-mark; Third algorithm: +- mark; 10 -7 : solid line. 



5.2. Choice of n. From the experiments in [11|,|12]. we have observed that ML(n)BiCGStab 
behaves more and more robust as n is increased. So, for an ill-conditioned problem, we would 
tend to suggest a large n for ML(n)BiCGStab. On the other hand, ML(n)BiCGStab minimizes 
||r-fc||2 once every n iterations. The convergence of a well-conditioned problem is usually accelerated 
by the minimization steps. So, when a problem is well-conditioned, we would suggest a small n. 

In [UCEO], it was suggested to fix s = 4 or 8 for the general use of IDR(s). This good idea also 
applies to ML(n)BiCGStab, namely, fixing n = 4 or 8 in its general use. 

We believe that the most powerfulness of ML(n)BiCGStab is in the solution of a sequence 
of linear systems. We once tested the first algorithm (see Algorithm 9.1 in |11| ) with n = 9 
and k = (see §5.31 for k) on the standard oil reservoir simulation test data called SPE9 and 
found that ML(n)BiCGStab reduced the total computational time by over 70% when compared 
to BiCGStab. A later test on SPE9 with Code #4 in p] showed that a 60% reduction in time 
can be reached. 

Code #4 is a design of automatic selection of the parameter n during the solution of a sequence 
of linear systems. Let tl and t2 denote the times to solve the previous and the current systems 
respectively. Then the basic idea behind Code #4 is: if tl > t2, then increase n to n + step when 
solving the next system; otherwise, decrease n to n — step. Here step is the search step size. 

We also plot the graphs of T conv (n) in Figure 15.21 to provide more information on how n affects 
the performance of ML(n)BiCGStab. 

5.3. Choice of uj. The standard choice for the w gri (k~\-i) in Algorithm l4. 1 1 (see Line 22) is w ffn (fe+i) = 
(Au fc ) H u fc /||Aufc||2. This choice of w Sn (fc+i) minimizes the 2-norm of r k = -uj gn ( k+1 )Au k + u k 
(Line 23), but sometimes can cause instability due to that it can be very small during an execution. 
The following remedy to guard 0J gn (k+l) aw ay from zero has been proposed in [5]: 

w g„(fe+l) = (Au fc ) H u fc /||Au fc |||; 
(5.1) p = (Au k ) H u k /(\\Au k \\ 2 ||u fc || 2 ); 

if \p\ < k, u g „(k+i) = KUJ 9n ( k+1 )/\p\; end 

where k is a user-defined parameter. See the numerical experiments in [Tllll] for more information 
about eqns (15.ip . 
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Figure 5.2. Graphs of 

T C onv( n ) against n. First algorithm: x-mark; Second 
algorithm: o-mark; Third algorithm: +-mark; BiCGStab: Solid line. BiCGStab 
took 2.07 and 10.64 seconds to converge for uim5940 and gc2534 respectively. 



6. Conclusions 

ML(n)BiCGStab is a powerful Krylov subspace method, especially in the solution of a sequence 
of linear systems with the parameter n dynamically chosen (see for detail) . This method has 
three algorithms. The first two can be found in and the third is new and is presented here 
as Algorithm 14.11 The third algorithm involves in its implementation and behaves as stable 
as the first algorithm, but converges faster than the first algorithm. Compared to the second 
algorithm, this third algorithm is more stable, but takes more time to converge. 



7. Appendix 

Algorithm 17.11 below is the preconditioned version of Algorithm 14.11 To avoid calling the index 
functions r n (k) and g n (k) every fc-iteration, we have split the fe-loop into a i-loop and a j-loop 
where k are related by (|2.ip with 1 < i < n, < j. 



Algorithm 7.1. ML(n)BiCGStabt with preconditioning 

1. Choose an initial guess xo and n vectors qi,q2, • • • , q n - 

2. Compute [fi, • • • ,f n -i] = M~ H A H [qi, • • • ,q n _i], r = b - Ax and g = r . 
Compute g = M _1 r , w = Ag , c = qf w , e = qf r . 



3. 


For j = 0,1,2,- •• 




4- 


For i = 1,2,--- ,n — 1 




5. 


Otjn+i — &jn+i—l / Cjn+i—1 j 




6. 


X-jn+i — "X-jn+i— 1 C^jn+iKjn+i- 


-1; 


7. 


Yjn+i — ^jn+i—l Otjn+iWjn+i- 


-1; 


8. 






9. 


If J > 1 




10. 


o(jn+i) j 
P(j-l)n+i ~ e jn+i/ %-l)n- 


H> 


11. 


n(jn+i) 

z w - r jn+i + ^ (j _ 1)n+i w (j _ 


l)n 


12. 


Sjn+i — P( J -_i) n+i g(j'-l)n+i/ 




13. 


For s = i + 1, ■ ■ ■ , n — 1 





Of oU n +i) _ _, , oijn+i) 
/0 P(j-l)n+i ~ U 'jP(j-l)n+i 



u 

15 

16. 
17. 

18 

19 

20 

21 

22 
23 

24 
25 
26. 
27. 

28 
29 
30 
31 
32 
33. 
34 
35. 
36. 
37. 
38. 
39 

40 
41 

42 

43. 

44 
45. 

46. 

47. 

48 
49 

50 

51 

52 
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o(jn+i) _ h _ /„ OX a(on+i) _ oijn+i) 

P(j-l)n+s ~ ^s+l £ w/ C(j-l) n+s , SO P(j_i) n+S — w iP(j-l) n +s 

z w — z w -t- p^.j)^ W(j_!) n+S ; 

gjn+i — gjn+i + P(j_i) n+s g(j-l)n+s/ 

End 

1 

gj'rt+i — z ui gjf n _)_j, 
LOj 

For s = 0, • • • , i — 1 

gjn+j = gjri+i + /3j n _|_ s gjn+s/ 

End 
Efee 

oO'n+i) _ _ffl_. / . 
Pjn ~~ A l l jn+i/^jn, 

Sjn+i Yjn+i ~l~ fijn ^>j n ' 

For s = 1, • • • , i — 1 

Pjn+a ^ = ~^s+lSjn+i/ 'cjn+s,' 

, n{jn+i) 

Sjn+i — gjn+i + Pj n _|_ s gjn+s; 

End 
End 

gj'ri+i — -M- gj'n+i; — Agjjj-j-j, 

End 

(Xjn+n — djn+n— l/Cjn+n— 1 > 

■^j'n+n — ^j'n+n— 1 "I - Q!jn+nKjn+n—l! 

Uj'n+n — *",j'rH-n— 1 Q-jn+fi " jn+n— 1 ; 

= (Auj„,+ n )^Uj„ + „/|| Auj„ + „|||; 

^-jn+n ^-jn+n ^j+l^jn+m 



o{jn+n) _ I ■ <¥ flO' n + n ) — _, , pU n + n ) 

P(j-l)n+n~ c jn+n/ L-(j-l)n+n, /(l P(j-l)n+„ " ^i+^-P {j-l)n+n 

Zto — Ijn+n "I" P^_i) n+n w Cj-l)n+n> 

_ o(jn+n) 
gjn+n — P( J _ 1 ) n+n g(i-l)n+n;' 

For s = 1, • • • , n — 1 

5(jn+n) _ w / _ qv 3(jn+n) _ , , o(jn+n) 

Pjn+s — q s +l z «)/ c in+s, /0 p s+ j n — LUj + ip g+ j n 

— _i_ 50' n + n ) 

. 5(j'n.+n) 

gjn+n — gj'n+n T" Pj n + S Kjn+si 

End 

Sjn+n — z u> Sjn+nj Sjn+n — -M- gj'n+n; 

- A - ^ - H 

End 



Matlab code of Algorithm 17.11 

1. function [x, err, iter, flag] = mlbicgstabt(A, x, b, Q, M, maxJt, tol, kappa) 
2. 
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3. % input: A: N-by-N matrix. M: N-by-N preconditioner matrix. 

4. % Q: N-by-n shadow matrix [qi ,•• • , q n ]. x: initial guess. 

5. % b: right hand side vector. maxJt: maximum number of iterations. 

6. % tol: error tolerance. 

7. % kappa: (real number) minimization step controller: 

8. % kappa = 0, standard minimization 

9. % kappa > 0, Sleijpen-van der Vorst minimization 

10. % output: x: solution computed, err: error norm, iter: number of iterations performed. 

11. % flag: = 0, solution found to tolerance 

12. % = 1, no convergence given max At iterations 

13. % = —1, breakdown. 

14. % storage: F: N x (n — 1) matrix. G,Q,W: N x n matrices. A, M: N x N matrices. 

15. % x,r, g_h, z,b: N x 1 matrices, c: 1 x n matrix. 
16. 

17. AT = size(A,2); n = size(Q,2); 

18. G = zeros(N,n); W = zeros(N,n); % initialize work spaces 

19. if n > 1, F = zeros(N, n — 1); end 

20. c = zeros (l,n); % end initialization 
21. 

22. iter = 0; //ag = 1; bnrm2 = norm(b); 

23. if 6nrm2 == 0.0, bnrml = 1.0; end 

24. r = 6 — A * x; err = norm(r) /bnrm2; 

25. if err < toZ, /Zag = 0; return, end 
26. 

27. if n > 1, F = M'\(A' * Q(:, 1 : n - 1)); end 

28. G(:, 1) = r; = M\r; 1) = A * gJi; c(l) = Q(:, 1)' * 1); 

29. if c(l) == 0, //a<7 = —1; return, end 

30. e = Q(:,l)'*r; 
31. 

32. for j = : max jit 

33. for i = 1 : n — 1 

34. alpha = e/c(i); x = x + alpha * gJi; r = r — alpha * W(:,i); 

35. err = norm(r)/bnrm2; iter = iter + 1; 

36. if err < to/, //ag = 0; return, end 

37. if iter >= maxJt, return, end 
38. 

39. e = Q(:,i + l)'*r; 

40. if j >= 1 

41. feeta = -e/c(i + 1); 

42. W(:,i + 1) = r + 6eta* W(:,i + 1); 

43. G(:,i + l) = beta*G(:,i + l); 

44. for s = i + 1 : n — 1 

45. feeta = -Q(:, s + 1)' * W(:, i + l)/c(s + 1); 

46. W(:, i + 1) = W(:, i + 1) + 6eta * W(:, s + 1); 

47. G(:, i + 1) = G(:, i + 1) + beta * G(:, s + 1); 

48. end 

49. G(:, i + 1) = W(:, i + 1) - G(:, i + 1)./ omega; 

50. for s = : i - 1 

51. beta = -F(:, s + 1)' * G(:, i + l)/c(s + 1); 
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52. G(:, i + 1) = G(:,i + 1) + beta * G(:, s + 1); 

53. end 

54. else 

55. beta = -F{:, 1)' * r/c(l); G(:, i + 1) = r + 6eto * G(:, 1); 

56. for s = 1 : i — 1 

57. 6eta = -F(:, s + 1)' * G(:, i + l)/c(s + 1); 

58. G(:, i + l) = G(:,i + 1) + beta * G(:, s + 1); 

59. end 

60. end 

61. g_h = M\G(:,i + l); W(:,i + 1) = A * gJi; 

62. c(i + l) = Q(:,i + l)'*W(:,i + l); 

63. if c(i + 1) == 0, flag = — 1; return, end 

64. end 

65. alpha = e/c(n); x = x + alpha * gJi; r = r — alpha * W(:,n); 

66. err = norm(r)/bnrm2; 

67. if err < tol, flag = 0; iter = iter + 1; return, end 

68. gJi = M\r; z = A* gJi; omega = z' * z; 

69. if omega == 0, flag = —1; return, end 

70. rho = z' * r; omega = rho/ omega; 

71. if kappa > 

72. rho = rhoj{norm(z) * norm(r)); absjom = abs(rho); 

73. if (abs_om < kappa) h (abs-om ~= 0) 

74. omega = omega* kappa/ abs-om; 

75. end 

76. end 

77. if omega == 0, flag = —1; return, end 

78. x = x + omega * gJi; r = r — omega * z; 

79. err = norm(r)/bnrm2; iter = iter + 1; 

80. if err < tol, flag = 0; return, end 

81. if iter >= maxJt, return, end 
82. 

83. e = Q(:, 1)' * r; beta = -e/c(l); 

84. W(:, 1) = r + feeta * W(:, 1); G(:, 1) = beta * G(:, 1); 

85. for s = 1 : n — 1 

86. 6eta = -Q(:, s + 1)' * W(:, l)/c(s + 1); 

87. W(:,l) = W(:,l) + beta*W(:,s + l); 

88. G(:,l) = G(:,1) + beta *G(:, s + 1); 

89. end 

90. G(:, 1) = W(:, 1) - G(:, 1). /omega; g_h = M\G(:, 1); 

91. W(:,l) = A*g_h; c{l) = Q {:,!)' *W {:,!); 

92. if c(l) == 0, //a^ = —1; return, end 

93. end 



A sample execution of ML(n)BiCGstabt 

1. N = 100; A = randn(N); M = randn(N); b = randn(N, 1); 

2. n = 10; kappa = 0.7; tol = 10~ 7 ; maxJt = 3 * N; 

3. Q = sign(randn(N, n)); x = zeros(N, 1); Q(:, 1) = 6 — ^4 * x; 

4. [x, err, iter, /Zag] = mlbicgstabt(A, x, b, Q, M, max-it, tol, kappa); 
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